Clustering and Classification of Satellite Imagery
Abstract
-RAQUEL # Unsupervised Clustering -RAQUEL any intro stuff about this topic, feell free to use references This section will use unsupervised clustering to split up given satellite images into several groups, with the goal of identifying meaningfull clusters that represent different geographic features such as water, farm land, forest or buildings.
## Loading required package: randomForest## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'randomForest'- Load a raster file, which has 4 layers
rasterJ<-brick("C:\\Users\\Owner\\Documents\\Portugal\\Sat_Image_Classification\\data\\data\\img\\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF")
#rasterJ<-brick("C:\\Users\\D059331\\Desktop\\DM GIC\\data\\img\\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF")
rasterJ## class : RasterBrick
## dimensions : 4238, 4239, 17964882, 4 (nrow, ncol, ncell, nlayers)
## resolution : 2.4, 2.4 (x, y)
## extent : -105092.3, -94918.74, -50088.84, -39917.64 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : C:\Users\Owner\Documents\Portugal\Sat_Image_Classification\data\data\img\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF
## names : J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.1, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.2, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.3, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.4
## min values : 0, 0, 0, 0
## max values : 255, 255, 255, 255- Plot map
plotRGB(rasterJ, 3,2,1)
- choose a sector for clustering
ext <- extent(-104294.4, -102964.5, -43623.48, -42742.44 )
sector <- crop(rasterJ, ext)
plotRGB(sector, 3, 2, 1)
- add normdiff layer for better vegetation recognition
sector[[5]] <- normdiff(sector)- plot each layer separately
plot(sector)
- Get summary statistics of layers. Note that there are no NAs
summary(sector)## J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.1
## Min. 23
## 1st Qu. 31
## Median 34
## 3rd Qu. 40
## Max. 176
## NA's 0
## J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.2
## Min. 26
## 1st Qu. 40
## Median 50
## 3rd Qu. 65
## Max. 255
## NA's 0
## J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.3
## Min. 8
## 1st Qu. 17
## Median 34
## 3rd Qu. 51
## Max. 194
## NA's 0
## J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.4 values
## Min. 4 -0.64885496
## 1st Qu. 10 -0.20000000
## Median 53 0.05645161
## 3rd Qu. 70 0.15384615
## Max. 191 0.76870748
## NA's 0 0.00000000- Input each of the five layers to the moving function. This creates windows of three by three pixels and calculates the mean and standard deviation of the pixels in these windows.Returns five rasters with these mean and standard deviation values filled in to the pixels.
means<-mov.fun(sector,3,mean, normalize=FALSE, verbose=TRUE)## >> Band 1 of 5 done.
## >> Band 2 of 5 done.
## >> Band 3 of 5 done.
## >> Band 4 of 5 done.
## >> Band 5 of 5 done.sds<-mov.fun(sector,3,sd, normalize=FALSE, verbose=TRUE)## >> Band 1 of 5 done.
## >> Band 2 of 5 done.
## >> Band 3 of 5 done.
## >> Band 4 of 5 done.
## >> Band 5 of 5 done.- convert output of above to brick and summarize to check for na values
stk<-stack(means,sds)
sectorMovings<-brick(stk)
summary(sectorMovings)## layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1 layer.1
## Min. 24.00000 27.88889 9.111111 4.555556 -0.51877395
## 1st Qu. 30.88889 40.44444 16.666667 10.111111 -0.21001302
## Median 34.55556 51.00000 35.555556 53.555556 0.05963125
## 3rd Qu. 40.22222 64.88889 50.888889 70.222222 0.16173571
## Max. 172.11111 255.00000 189.222222 171.111111 0.73729629
## NA's 1840.00000 1840.00000 1840.000000 1840.000000 1840.00000000
## layer.1.1.2 layer.2.1.2 layer.1.2.2 layer.2.2.2 layer.2
## Min. 0.0000000 0.000000 0.000000 0.000000 0.000000e+00
## 1st Qu. 0.9718253 1.666667 1.322876 1.414214 2.197058e-02
## Median 1.4142136 2.877113 2.438123 2.915476 4.184848e-02
## 3rd Qu. 2.8037673 6.009252 5.967505 7.270565 6.507355e-02
## Max. 53.2846132 85.835275 64.031242 60.055761 3.334456e-01
## NA's 1840.0000000 1840.000000 1840.000000 1840.000000 1.840000e+03-We see that there are equal NA values in each layer We get dimension information to determine number of cells in total in each layer
sectorMovings## class : RasterBrick
## dimensions : 367, 555, 203685, 10 (nrow, ncol, ncell, nlayers)
## resolution : 2.4, 2.4 (x, y)
## extent : -104295.5, -102963.5, -43623.24, -42742.44 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1, layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2, layer.2
## min values : 24.0000000, 27.8888889, 9.1111111, 4.5555556, -0.5187739, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000
## max values : 172.1111111, 255.0000000, 189.2222222, 171.1111111, 0.7372963, 53.2846132, 85.8352751, 64.0312424, 60.0557611, 0.3334456-We see that NA values are a very small percentage of total cells (less than one percent) of all cells#they are not visible in the plots
sectorMovings## class : RasterBrick
## dimensions : 367, 555, 203685, 10 (nrow, ncol, ncell, nlayers)
## resolution : 2.4, 2.4 (x, y)
## extent : -104295.5, -102963.5, -43623.24, -42742.44 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1, layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2, layer.2
## min values : 24.0000000, 27.8888889, 9.1111111, 4.5555556, -0.5187739, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000
## max values : 172.1111111, 255.0000000, 189.2222222, 171.1111111, 0.7372963, 53.2846132, 85.8352751, 64.0312424, 60.0557611, 0.3334456plot(sectorMovings)
-use trim to remove NA values -from documentation: Trim (shrink) a Raster* object by removing outer rows and columns that all have the same value (e.g. NA). -Check output for NA values
sectorMovingsTrim<-trim(sectorMovings, padding=0, values=NA)-No NA values in any layers -Plote output of trim to see if any issues
sectorMovingsTrim<-trim(sectorMovings, padding=0, values=NA)
summary(sectorMovingsTrim)## layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1 layer.1
## Min. 24.00000 27.88889 9.111111 4.555556 -0.51877395
## 1st Qu. 30.88889 40.44444 16.666667 10.111111 -0.21001302
## Median 34.55556 51.00000 35.555556 53.555556 0.05963125
## 3rd Qu. 40.22222 64.88889 50.888889 70.222222 0.16173571
## Max. 172.11111 255.00000 189.222222 171.111111 0.73729629
## NA's 0.00000 0.00000 0.000000 0.000000 0.00000000
## layer.1.1.2 layer.2.1.2 layer.1.2.2 layer.2.2.2 layer.2
## Min. 0.0000000 0.000000 0.000000 0.000000 0.00000000
## 1st Qu. 0.9718253 1.666667 1.322876 1.414214 0.02197058
## Median 1.4142136 2.877113 2.438123 2.915476 0.04184848
## 3rd Qu. 2.8037673 6.009252 5.967505 7.270565 0.06507355
## Max. 53.2846132 85.835275 64.031242 60.055761 0.33344559
## NA's 0.0000000 0.000000 0.000000 0.000000 0.00000000plot(sectorMovingsTrim)
-It looks good. The moving function may have created NA values around edges, maybe where a full three by three grid could not be made.
-determine kmeans seed values
-Normalize the ten layer stack because the range of values and variation differ between the RGB and infrared bands the NDVI. The five layers measuring the standard deviation also have a different range and variation level.
sectorMovingsTrimNorm<-raster.scale(sectorMovingsTrim)
sectorMovingsTrimNorm## class : RasterBrick
## dimensions : 365, 553, 201845, 10 (nrow, ncol, ncell, nlayers)
## resolution : 2.4, 2.4 (x, y)
## extent : -104293.1, -102965.9, -43620.84, -42744.84 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## data source : in memory
## names : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1, layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2, layer.2
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1-run kmeans on the normalized ten layer stack.
performKMeans(sectorMovingsTrim,6)## Warning: Quick-TRANSfer stage steps exceeded maximum (= 10092250)
-Assess the clusters we have found -RAQUEL any writing about how we assess clusters (generally, or in R)
-Add the the supervised report here